Numerical investigation of the influence of the history on the local structure of glasses 
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By means of molecular dynamics simulations and the Vorono'f tessellation, we study the influence of 
the history on the low temperature characteristics of soft sphere and silica glasses. The quench from 
the liquid is interrupted at an intermediate temperature Ti for a given relaxation time, and then 
the cooling process is continued down to OK. The local structure at OK depends on the temperature 
Ti and the effect is larger for T close to the glass transition temperature T g . This dependence, 
coherent with recent results, is expected in the strong glass former where the characteristics of a 
particular state depend on its history. In the soft-sphere case, because of crystallization effects, the 
dependence of the local structure of the glassy samples on their history can only be detected in the 
supercooled liquid region. 
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I. INTRODUCTION 

Experimentally, from the study of slow dynamics in 
disordered solids, it is now well established that the time 
evolution of certain systems depends on the manner these 
systems were prepared. This phenomenon called usually 
"aging" has been observed in structural (polymer) glasses 
0, spin-glasses 0] or orientational glasses S. Recently 
theoretical work concerning aging effects in a structural 
Lennard- Jones glass has been published J3J. Through 
classical molecular dynamics (MD) simulations the au- 
thors show that the dynamical behavior of this system 
depends on the initial temperature before the quench 
and exhibits the same scaling features than the one ob- 
served in spin-glasses. Several theoretical explanations 
have been proposed to account for this apparent univer- 
sal behavior of the out-of-equilibrium dynamics in various 
glasses ||. 

Since it has been shown that aging effects occur in a 
model glass, we want here to address this question from 
the structural point of view by using a combination of 
classical microcanonical MD calculations and the Voronoi 
tessellation. An in-depth description of our Voronoi tes- 
sellation scheme can be found in a previous paper Q . For 
the sake of generality, we consider two very different sys- 
tems, a soft sphere glass and a silica glass, representative 
of amorphous metals and vitreous oxides, respectively. 



II. METHODOLOGY 

The soft sphere system is a very simple monoatomic 
model (without any particular experimental counterpart) 
described by the purely repulsive potential introduced by 
Laird and Schober (LS) Ij. This potential is basically a 
modified inverse 6 th power potential similar to those very 
often used in computer simulations of liquids and glasses. 
Here we use the same parameters (arbitrarily equal to 



the Lennard- Jones parameters of Argon) as in previous 
studies HH and consequently the glass transition tem- 
perature T g and melting temperature T m are respectively 
T g ~ 10.3K § and T m ~ 22. 7K §. In the present MD 
calculations we have considered, as before |^^|, 1000 soft 
spheres in a cubic box of edge length 34.05A. The sil- 
ica system is described by the more sophisticated po- 
tential hrst introduced by van Beest, Kramer and van 
Santen (BKS) M] and justified by ab initio calculations. 
Though designed originally for the crystalline phases of 
silica, it has been shown to describe quite well the struc- 
tural [jllj and vibrational |l2] properties of vitreous sil- 
ica. In that case also we make use of the same parameters 
as in a previous study |13| in which we have performed 
microcanonical MD calculations (treating the Coulomb 
interactions using the Ewald summation method pi). 
The Voronoi tessellation procedure has been straightfor- 
wardly adapted to this two component system and it is 
worth noting that we did not include any weighing factor 
to distinguish silicon from oxygen atoms |L3| |. The esti- 
mated value for the glass transition temperature (ll] [H| 
is T g ~ 3500K, subsequently larger than the experimen- 
tal one partly because of the very fast cooling rate (2.3 
10 14 K/s) imposed by numerical calculation limitations. 
But the value of T g is consistent with its dependence with 
the cooling rate proposed by Vollmayr et al. [fl"l| in which 
a Vogel-Fulcher dependence of the relaxation time r* of 
the system on the temperature is assumed. As in [fl3"f , our 
silica system consists of 216 silicon and 432 oxygen atoms 
confined in a cubic box of edge 2 1.48 A corresponding to 
a density close to the experimental one. In the following 
the soft sphere and silica systems will be called LS and 
BKS, respectively. 

The purpose of the present work is to study the influ- 
ence of the sample history on the low temperature struc- 
tural properties of our two model glasses. To achieve 
this goal we follow in both cases a very simple modus 
operandi, which is illustrated in Fig. 1 in the LS case. 
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We start from well equilibrated liquid samples at a well 
denned temperature Ti (point A). We have chosen Te = 
45K in the LS case and 7000K in the BKS case. At this 
point two paths can be followed : either the liquid is 
quenched to an intermediate temperature T (point B) 
then relaxed in a microcanonical ensemble during a time 
t (point C) and finally quenched to zero temperature 
(point D) or the liquid is quenched to OK directly (point 
C) relaxed during r and finally cooled to OK (point D). 
This last cooling procedure is necessary because of the 
slight temperature increase following the direct quench 
to OK since our system evolves at constant energy. From 
a structural point of view we have checked that this short 
quench is not really necessary but we have included it 
anyway in the second path to have a unique final state 
(point D). 

III. SOFT-SPHERE GLASS 

In the LS case, the quenching rate 7 has been chosen 
equal to the "ideal" quenching rate of 10 12 K/s for which 
the glass stability at OK compared to the crystalline state 
has been shown to be maximum [|| and the local struc- 
ture at point D has been studied using the fraction of 
pentagonal faces f$ of the Vorono'i cell attached to each 
particle. We have previously shown that this parameter 
is a very sensitive tool : a high value of f$ (typically 
0.45) is a sign of strong icosahedral order, characteristic 
of amorphous metallic glasses, while a small value of /g 
(<0.2) indicates the onset of crystal nucleation. Indeed, 
depending on t, one has to be careful in order to dis- 
tinguish between pure glassy state relaxation and recrys- 
tallization effects. Here we have considered two different 
values of r : 100 and 300 ps. These values of r should 
be compared to the value of the typical relaxation time 
of the system, t* . This relaxation time depends on the 
quenching rate. It is generally admitted that the sys- 
tem falls out of equilibrium when its relaxation time is 
on the order of the time scale of the cooling process i.e. 
t* as Tgj^j as 32 ps where 7 is the effective cooling rate 
deduced from Fig.l. The values of r used in our study 
are therefore large enough to permit a complete relax- 
ation of the systems for Tj > T g . Of course this is no 
more the case in the solid phase. 

Concerning T,-, for each sample we saved configurations 
(positions and velocities) along the quench at fixed time 
intervals corresponding approximately to a difference of 
temperature of 2K. The corresponding results have then 
been averaged over ten different samples generated from 
ten independent liquid configurations. It is interesting 
to note in Fig.l that the amplitude of the temperature 
fluctuations during the relaxation procedure follows the 
usual VT behavior and are almost inexistent at very low 
temperature. 

In Fig. 2 we show the variation of /g as a function of T for 



t = 100 and r = 300 ps. The error bars show the tem- 
perature dependent dispersion of our statistical ensemble 
containing ten samples. The first information concerns 
the high temperature behavior. Fig. 2 points out that 
if Ti is above a temperature T sup as I9K (lower than 
T m ), the structure obtained at OK is the same whichever 
path the system followed during the quench. This struc- 
ture is symbolically represented in Fig. 2 by the horizontal 
dashed line. Indeed the points for T$ > I9K lie nicely on 
this line as well as the extrapolated value for T = 
(due to averaging, the value Tj = can not be reached 
exactly). This structure is without any doubt a glassy ar- 
rangement and exhibits a value of /g close to 0.45 which 
is a sign of strong icosahedral local order. This first in- 
formation tells us that the low temperature structure fol- 
lowing the quench will be the same whatever the initial 
liquid temperature is, as long as it is above T sup . In the 
picture representing the whole system as a point evolv- 
ing on a complex multi-valleyed energy "surface" [fl5"|| , it 
means that as long as the temperature Ti is above T sup , 
the kinetic energy of the system is greater than all the en- 
ergy barriers on that surface. Therefore if Ti > T sup the 
system explores the whole configuration space and the 
low temperature sample does not depend on T. From a 
purely practical point of view it means also that quenches 
from the liquid state can be started at temperatures just 
above the melting temperature without the system re- 
membering the characteristics of this liquid state. 
The curves in Fig. 2 depart from the horizontal dashed 
line (which represents in a sense the "reference" struc- 
ture) for T below as 19K. It is interesting to note that 
the two curves corresponding to the two values of r are 
superimposed down to as 15K. Then these two curves 
split and exhibit a minimum for T as 10K. This mini- 
mum is close to 0.38 for r = lOOps and is much lower for 
r = 300ps ( /g as 0.12). This is a direct consequence of 
the relaxation, the system after lOOps being on its way to 
the equilibrium configuration reached after 300ps. What 
is the meaning of this minimum ? As said earlier a small 
value of /g is a sign of the onset of crystallization. More- 
over looking at the results obtained for the individual 
samples it appears that some systems exhibit values of 
/g as small as I0 -2 indicating an almost perfect crys- 
talline character. This can also be verified by direct vi- 
sualization of these samples. These results indicate that 
if T is chosen close to the glass transition temperature T g 
{T g = 10. 3K), the propensity to crystallize is maximum. 
This means that our model glass exhibits a maximum 
of instability with respect to the crystalline state when 
its temperature is close to T g . Below T g , /g increases 
again and the low temperature glass stability is recov- 
ered as shown previously for this quenching rate Q . The 
fact that the tendency to crystallize increases when T 
increases up to T g can be explained using the same kinds 
of arguments than the ones developed in js| . When T in- 
creases, the local density fluctuations increase and there- 
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fore the probability to find crystal germs in the samples 
increases. Of course, as in [B, this explanation remains 
qualitative since it does not take into account the ability 
of the existent germs to grow. To explain the existence 
of a minimum around T g in the /s curve i.e. that the 
stability increases when Tj increases above T g , one can 
argue that the crystal germ nucleation is now outweighed 
by the fact that the energy difference between the super- 
cooled state and the crystal state increases, making the 
crystallization process more difficult to happen. Let us 
discuss again in terms of relaxation times. When crys- 
tal nucleation occurs the nucleation rate has to be taken 
into account introducing an extra time scale, r„(T). r„ is 
a very rapidly varying, non monotonous function of the 
temperature: it first decreases as T is lowered from the 
liquid state, due to the decrease in the free energy barrier 
and then increases again due to the liquid becoming more 
sluggish. This variation is shown in the so-called "TTT" 
diagrams The curves in Fig. 2 can then be inter- 

preted in the following way: the waiting time t = lOOps 
is always smaller than r n while r = 300ps becomes larger 
than r n between 12.5 and 15K. Hence the samples relaxed 
at Ti < Ti„f 13K exhibit crystal nucleation. The max- 
imum crystalline order is reached when the nucleation 
rate is maximum i.e. around T g |16|. 

When the curves in Fig. 2 depart from the horizon- 
tal dashed line, it means that the glass "remembers" its 
past history and has a different structure depending on 
the path (ABC'D or ACD) which it followed during the 
quench. Nevertheless most of these differences are due 
to recrystallization effects which are interesting but are 
not aging effects (the term aging was in fact introduced 
to "distinguish glassy state relaxation from other time- 
dependent processes such as recrystallization" To 
make this distinction we report in Fig. 3 the variation for 
r = 300ps of /s as a function of Ti (full circles) but also 
as a function of Tf (open squares). Tf is the temperature 
of the system at point C in Fig.l. It corresponds to the 
final temperature after the waiting time r and can be 
compared directly to the temperature used in the work 
done on the Lennard- Jones glass As can be seen 
in Fig. 3, these temperatures are almost identical which 
permits to compare the /5S in the system at C and D. 
The difference between these two states is an additional 
quench to lower the temperature from Tf to OK. In fact 
the time spent during this quench can be viewed as an 
extra relaxation time. Therefore if the system evolves 
towards a crystalline state this supplemental relaxation 
time should lead to a decrease of the /5S at D compared 
to the results obtained at C. In Fig. 3 this is the case for 
temperatures between 2K and the crossover temperature 
Tinf. This indicates that the departure from the refer- 
ence structure in that temperature range is simply due 
to crystallization effects. In the liquid region for temper- 
atures above T sup the smaller values of /s at Tf are a 
consequence of thermal effects 071 and therefore it is ex- 



pected that once the system is cooled down the values of 
/5 increase towards the values obtained in the reference 
structure. Let us now examine the last temperature in- 
terval which goes from Ti n f to T sup . In Fig. 3 we see that 
the additional relaxation time induces an increase of the 
/5S in that temperature region (as an example at 15K, 
f 5 is equal to 0.28 at C while at D f 5 = 0.42). There- 
fore no onset of crystallization is present in the samples. 
On the contrary the additional quench makes them even 
more amorphous. Nevertheless in that temperature in- 
terval the curve of /s as a function of Ti is not super- 
imposed with the horizontal dashed line. Therefore the 
glassy structures obtained are different from the reference 
structure indicating the effect of the history on the local 
structure of the amorphous phase. The differences are 
admittedly small but they exist and are larger than the 
error bars in that temperature range. This effect exists 
between Ti„/ and T sup i.e. roughly in the supercooled liq- 
uid region. It may exist also at lower temperatures but to 
us it will remain hidden by the recrystallization effects. 
In that temperature region the system does not explore 
the whole configuration space anymore and therefore the 
glassy systems obtained are different from the ones ob- 
tained at higher temperatures. Locally they appear less 
disordered (^5 is smaller than the value obtained for the 
reference structure) because the system has not access to 
the full variety of possible configurations but it remains 
amorphous. In fact we have a picture of the configura- 
tional space corresponding to Ti in which the crystalline 
ground state is not accessible and which leads to glassy 
systems structurally different from the reference struc- 
ture. This difference is quantified in Fig. 3. 



IV. SILICA GLASS 

In the BKS case, the quenching rate 7 has been cho- 
sen equal to 2.3 10 14 K/s and a unique relaxation time 
r = 42ps (corresponding to 60000 timesteps) has been 
considered. This choice for r results from a balance be- 
tween an upper limit imposed by our computing facilities 
and a lower limit necessary to obtain detectable effects, 
even after averaging over ten independent samples, as we 
did here. Also this value of r is nearly three times larger 
than the typical relaxation time of the system, r* 15ps, 
which permits a complete relaxation of the samples dur- 
ing the waiting time for Tj > T g . The difficulty in obtain- 
ing these results is rewarded by the certitude that in the 
BKS case no crystallization tendency has ever been seen 
Pi O . To detect the structural modifications in the 
case of vitreous silica we can no more make use of the /s 
parameter since it has a significance only for amorphous 
monoatomic glasses with spherical symmetry. Here we 
have chosen to present the standard deviations o~o and 
o~Si of the Vorono'i cell volumes for the oxygen and silicon 
atoms over the whole sample. In fact we have calculated 
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the dimensionless quantities So = °o/ < Vo > and 
Ssi = a si/ < Vsi > in which the standard deviation 
has been divided by the corresponding averaged cell vol- 
ume. The standard deviation of the cell volumes is a 
direct measure of the local density fluctuations around a 
given atom. As shown in a recent study this quan- 
tity decreases when lowering the temperature from the 
liquid state and it is possible to extrapolate it to zero 
at OK as if the system would try to reach a crystalline 
state. This decrease slows down when T becomes lower 
than the glass transition temperature T g , and finally the 
standard deviations saturate to a non zero value at OK, 
typical of glassy disorder. Therefore the lower the cr's 
are, the closer to an "ideal" glassy state the considered 
system is (as long as there is no crystallization). The 
variation of S as a function of Ti is plotted in Fig. 4. 
Note that Egi is smaller than So which reflects a larger 
disorder around the oxygen atoms, in agreement with the 
fact that the tetrahedral arrangement around the silicon 
atoms is quite well preserved, even in the liquid or amor- 
phous states. In the case of silica also, the data obtained 
for a sufficiently large temperature (Tj > T sup w 6000K) 
correspond to the one obtained for = OK (represented 
by the horizontal dashed lines in the figure) within the er- 
ror bars. The small but significant data departure from 
the dashed line in the whole range < Tj < T sup is 
maximal when Ti is chosen close to T g , similar to the 
behavior observed in the soft-sphere glass. How can this 
departure be interpreted in the case of silica in which no 
crystallization occurs ? For a strong glass former, at suf- 
ficiently low temperature, the characteristics of a sample 
depend on the history of the sample (quenching rate, re- 
laxation time) as shown previously plj j . Therefore we 
expect that the local structure of our silica glass depends 
on Ti and this is indeed the case as shown quantitatively 
in Fig. 4. Nevertheless in that point of view the shape 
of the curve needs to be explained: why is there a min- 
imum close to T g ? For temperatures below T g , in the 
solid phase, the structural relaxation time exceeds by far 
40ps and therefore our waiting time r is not sufficient 
to permit a complete relaxation of the systems and the 
situation gets worse as the temperature decreases. This 
explains why the curves in Fig. 4 converge towards the 
reference lines at low temperature. At T g , r* « 15ps 
and the waiting time is sufficient to let the systems re- 
lax closer the underlying ideal glassy structure (with a 
smaller S). In terms of energy landscapes, at T g the sys- 
tem has enough kinetic energy to overcome the barriers 
separating the metastable state characteristic of the con- 
sidered quenching rate from the lower lying states, during 
the waiting time. At a higher temperature, the system 
explores more phase space during the time r and there- 
fore gets trapped after the quench in a state close to the 
metastable state characteristic of the cooling rate. With 
still higher temperatures, the system explores the same 
energy landscape during the waiting time and the low 



temperature structure (after the quench) is just a reflec- 
tion of the fast cooling rate. The minimum in the curves 
in Fig. 4 indicates that the structure obtained by stop- 
ping the quench to perform a long relaxation period at 
a temperature close to T g is closer to the "ideal" glassy 
structure, since it corresponds to lower values for the 
standard deviations of the cell volumes than in the ref- 
erence system (quenched at the same speed, but without 
interruption). This is consistent with the fact that, for a 
slower quench, the glass transition temperature would be 
smaller and so would be the standard deviations of the 
cell volumes (which have been shown to decrease with 
decreasing temperature in the liquid phase £31). In prac- 
tice, to obtain the same result, it is however more efficient 
to stop the quench and perform a relaxation at T g , which 
is the procedure followed by the authors in than to 
perform a slower purely linear quench. 

V. CONCLUSION 

This study combining MD simulations and the Voronoi 
tessellation in model glasses provides an other point of 
view in the study of the influence of the glass history 
on its properties. Here we tackle this question from the 
structural point of view. First of all the study as a func- 
tion of temperature of the local structure in soft-sphere 
samples has shown that this model glass exhibits a max- 
imum instability with respect to the crystalline state at 

We have then shown that the structure of this model 
glass depends on its history. Indeed when a system is 
cooled to an intermediate temperature, relaxed and then 
quenched to zero temperature (and if it stays amorphous) 
its structure is different from the structure obtained by a 
direct quench to OK followed by the same relaxation pe- 
riod if the intermediate temperature is chosen roughly in 
the supercooled liquid region. Secondly we did the same 
kind of study in a model silica glass where the crystal- 
lization phenomenon is absent. In that case we find also, 
as expected, a dependence of the local structure on the 
glass history. This dependence is maximum when the in- 
termediate temperature is close to T g . 
These results show that the Voronoi cell characteristics 
permit to clearly distinguish recrystallization from "true" 
glassy state relaxations in the case of the soft-sphere 
glass. In the case of silica an extension of this work would 
be to check if the change of the local structure observed 
close to T g is detectable in other quantities. 
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FIG. 1. Example of the quenching procedure used in this 
study: a relaxed liquid (point A) is quenched to OK following 
either the path ACD or the path ABC'D. 

Figure 2 
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FIG. 4. Variation of the reduced standard deviations 
E = a I < V >, for the oxygen (o) and the silicon (•) atoms 
in our silica system, as a function of the intermediate tem- 
perature Ti for r = 42 ps. The horizontal lines show the 
corresponding values in the reference structure at OK and are 
guides for the eye. 



FIG. 2. Variation of the fraction of pentagonal faces, /s, 
as a function of the intermediate temperature T» for r = 100 
ps ( o) and r = 300ps (•). The horizontal dashed line is a 
guide for the eye and indicates the value of /s in the reference 
structure at OK. 

Figure 3 
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FIG. 3. f 5 as a function of Ti (•) and T f (□). T s is the final 
temperature after a relaxation time r = 300ps (temperature 
at point C in Fig. 1). The horizontal dashed line is a guide for 
the eye and indicates the value of /s in the reference structure 
at OK. 
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